home *** CD-ROM | disk | FTP | other *** search
/ SGI Hot Mix 17 / Hot Mix 17.iso / HM17_SGI / research / lib / int_2d.pro < prev    next >
Text File  |  1997-07-08  |  13KB  |  272 lines
  1. ; $Id: int_2d.pro,v 1.8 1997/01/15 03:11:50 ali Exp $
  2. ;
  3. ; Copyright (c) 1995-1997, Research Systems, Inc.  All rights reserved.
  4. ;       Unauthorized reproduction prohibited.
  5. ;+
  6. ; NAME:
  7. ;       INT_2D
  8. ;
  9. ; PURPOSE:
  10. ;       This function computes the double integral of a bivariate
  11. ;       function F(x,y) using either a dy-dx or dx-dy order of integration.
  12. ;
  13. ; CATEGORY:
  14. ;       Numerical Analysis.
  15. ;
  16. ; CALLING SEQUENCE:
  17. ;       Result = INT_2D(Fxy, AB_Limits, PQ_Limits, Pts)
  18. ;
  19. ; INPUTS:
  20. ;       Fxy:  A scalar string specifying the name of a user-supplied
  21. ;             IDL function that defines the bivariate function to be
  22. ;             integrated. The function must accept x & y and return
  23. ;             a scalar result.
  24. ;       AB_Limits:  A two-element vector containing the lower, A, and
  25. ;                   upper, B, limits of integration. These limits correspond
  26. ;                   to x when using a dy-dx order of integration or they
  27. ;                   correspond to y when using a dx-dy order of integration. 
  28. ;       PQ_Limits:  A scalar string specifying the name of a user-supplied 
  29. ;                   IDL function that defines the lower, P, and upper, Q, 
  30. ;                   limits of integration. These limits correspond to y when
  31. ;                   using a dy-dx order of integration or they correspond to
  32. ;                   x when using a dx-dy order of integration. They must 
  33. ;                   accept x (for a dy-dx order of integration) or y (for a
  34. ;                   dx-dy order of integration) and return a two-element
  35. ;                   vector result.
  36. ;       Pts:  The number of transformation points used in the
  37. ;             computation. Possible values are: 6, 10, 20, 48 or 96.
  38. ;
  39. ; KEYWORD PARAMETERS:
  40. ;       DOUBLE: If set to a non-zero value, computations are done in
  41. ;               double precision arithmetic.
  42. ;       ORDER:  A scalar value of either 0 or 1. If set to 0, the integral 
  43. ;               is computed using a dy-dx order of integration. If set to 1,
  44. ;               the integral is computed using a dx-dy order of integration.
  45. ;
  46. ; EXAMPLE:
  47. ;       Compute the double integral of the bivariate function
  48. ;       F(x, y) = y * cos(x^5) over the region: [0.0, x^2] for y and 
  49. ;                                               [0.0, 2.0] for x using 
  50. ;                                               a dy-dx order of integration.
  51. ;
  52. ;       ;Define the limits of integration for y as a function of x.
  53. ;         function PQ_Limits, x
  54. ;           return, [0.0, x^2]
  55. ;         end
  56. ;
  57. ;       ;Define the limits of integration for x.
  58. ;         AB_Limits = [0.0, 2.0]
  59. ;
  60. ;       ;Integrate with 48 and 96 point formulas using a dy-dx order of 
  61. ;       ;integration and double-precision arithmetic.
  62. ;         print, INT_2D('Fxy', AB_Limits, 'PQ_Limits', 48, /double) 
  63. ;         print, INT_2D('Fxy', AB_Limits, 'PQ_Limits', 96, /double)
  64. ;
  65. ;       INT_2D with 48 transformation points yields:    0.055142668
  66. ;       INT_2D with 96 transformation points yields:    0.055142668
  67. ;
  68. ;       Compute the double integral of the bivariate function
  69. ;       F(x, y) = y * cos(x^5) over the equivalent region using a dx-dy order
  70. ;       of integration.
  71. ;
  72. ;       ;Define the limits of integration for x as a function of y.
  73. ;         function PQ_Limits, y
  74. ;           return, [sqrt(y), 2.0]
  75. ;         end
  76. ;
  77. ;       ;Define the limits of integration for y.
  78. ;         AB_Limits = [0.0, 4.0]
  79. ;
  80. ;       ;Integrate with 48 and 96 point formulas using a dx-dy order of
  81. ;       ;integration and double-precision arithmetic.
  82. ;         print, INT_2D('Fxy', AB_Limits, 'PQ_Limits', 48, /double, /order)
  83. ;         print, INT_2D('Fxy', AB_Limits, 'PQ_Limits', 96, /double, /order)
  84. ;
  85. ;       INT_2D with 48 transformation points yields:    0.055142678
  86. ;       INT_2D with 96 transformation points yields:    0.055142668
  87. ;
  88. ;       The exact solution (9 decimal accuracy) yields: 0.055142668
  89. ;
  90. ; PROCEDURE:
  91. ;       INT_2D.PRO computes the double integral of a bivariate function
  92. ;       using iterated Gaussian Quadrature. The algorithm's transformation
  93. ;       data is provided in tabulated form with 15 decimal accuracy.
  94. ;
  95. ; REFERENCE:
  96. ;       Handbook of Mathematical Functions
  97. ;       U.S. Department of Commerce
  98. ;       Applied Mathematics Series 55
  99. ;
  100. ; MODIFICATION HISTORY:
  101. ;       Written by:  GGS, RSI, January 1994
  102. ;       Modified:    GGS, RSI, September 1994
  103. ;                    Added 96 point transformation data.               
  104. ;                    Added DOUBLE keyword. Replaced nested FOR loop with
  105. ;                    vector operations resulting in faster execution.
  106. ;       Modified:    GGS, RSI, June 1995
  107. ;                    Added ORDER keyword allowing a dx-dy order of integration.
  108. ;       Modified:    GGS, RSI, April 1996
  109. ;                    Modified keyword checking and use of double precision.
  110. ;-          
  111.  
  112. FUNCTION Int_2D, Fxy, AB_Limits, PQ_Limits, Pts, Double = Double, Order = Order
  113.  
  114.   ON_ERROR, 2
  115.  
  116.   if N_ELEMENTS(AB_Limits) ne 2 then $
  117.     MESSAGE, "AB_Limits parameter must be a two-element vector."
  118.  
  119. ;Tabulated transformation data with 15 decimal accuracy.
  120. if Pts eq 6 then begin
  121.   Ri    = DBLARR(Pts)          &  Wi    = DBLARR(Pts)
  122.   Ri[0] = 0.932469514203152d   &  Wi[0] = 0.171324492379170d
  123.   Ri[1] = 0.661209386466265d   &  Wi[1] = 0.360761573048139d
  124.   Ri[2] = 0.238619186083197d   &  Wi[2] = 0.467913934572691d
  125.   Ri[INDGEN(Pts/2) + (Pts/2)]  = -Ri[(Pts/2) - INDGEN(Pts/2) -1]
  126.   Wi[INDGEN(Pts/2) + (Pts/2)]  =  Wi[(Pts/2) - INDGEN(Pts/2) -1]  
  127. endif else if Pts eq 10 then begin
  128.   Ri    = DBLARR(Pts)          &  Wi    = DBLARR(Pts)
  129.   Ri[0] = 0.973906528517172d   &  Wi[0] = 0.066671344308688d
  130.   Ri[1] = 0.865063366688985d   &  Wi[1] = 0.149451349150581d
  131.   Ri[2] = 0.679409568299024d   &  Wi[2] = 0.219086362515982d
  132.   Ri[3] = 0.433395394129247d   &  Wi[3] = 0.269266719309996d
  133.   Ri[4] = 0.148874338981631d   &  Wi[4] = 0.295524224714753d
  134.   Ri[INDGEN(Pts/2) + (Pts/2)]  = -Ri[(Pts/2) - INDGEN(Pts/2) -1]
  135.   Wi[INDGEN(Pts/2) + (Pts/2)]  =  Wi[(Pts/2) - INDGEN(Pts/2) -1]
  136. endif else if Pts eq 20 then begin
  137.   Ri     = DBLARR(Pts)         &  Wi     = DBLARR(Pts)
  138.   Ri[0]  = 0.993128599185094d  &  Wi[0]  = 0.017614007139152d
  139.   Ri[1]  = 0.963971927277913d  &  Wi[1]  = 0.040601429800386d
  140.   Ri[2]  = 0.912234428251325d  &  Wi[2]  = 0.062672048334109d
  141.   Ri[3]  = 0.839116971822218d  &  Wi[3]  = 0.083276741576704d
  142.   Ri[4]  = 0.746331906460150d  &  Wi[4]  = 0.101930119817240d 
  143.   Ri[5]  = 0.636053680726515d  &  Wi[5]  = 0.118194531961518d
  144.   Ri[6]  = 0.510867001950827d  &  Wi[6]  = 0.131688638449176d
  145.   Ri[7]  = 0.373706088715419d  &  Wi[7]  = 0.142096109318382d
  146.   Ri[8]  = 0.227785851141645d  &  Wi[8]  = 0.149172986472603d
  147.   Ri[9]  = 0.076526521133497d  &  Wi[9]  = 0.152753387130725d
  148.   Ri[INDGEN(Pts/2) + (Pts/2)]  = -Ri[(Pts/2) - INDGEN(Pts/2) -1]
  149.   Wi[INDGEN(Pts/2) + (Pts/2)]  =  Wi[(Pts/2) - INDGEN(Pts/2) -1]
  150. endif else if Pts eq 48 then begin
  151.   Ri     = DBLARR(Pts)         &  Wi     = DBLARR(Pts)
  152.   Ri[0]  = 0.998771007252426d  &  Wi[0]  = 0.003153346052305d
  153.   Ri[1]  = 0.993530172266350d  &  Wi[1]  = 0.007327553901276d
  154.   Ri[2]  = 0.984124583722826d  &  Wi[2]  = 0.011477234579234d
  155.   Ri[3]  = 0.970591592546247d  &  Wi[3]  = 0.015579315722943d
  156.   Ri[4]  = 0.952987703160430d  &  Wi[4]  = 0.019616160457355d
  157.   Ri[5]  = 0.931386690706554d  &  Wi[5]  = 0.023570760839324d
  158.   Ri[6]  = 0.905879136715569d  &  Wi[6]  = 0.027426509708356d
  159.   Ri[7]  = 0.876572020274247d  &  Wi[7]  = 0.031167227832798d
  160.   Ri[8]  = 0.843588261624393d  &  Wi[8]  = 0.034777222564770d
  161.   Ri[9]  = 0.807066204029442d  &  Wi[9]  = 0.038241351065830d
  162.   Ri[10] = 0.767159032515740d  &  Wi[10] = 0.041545082943464d
  163.   Ri[11] = 0.724034130923814d  &  Wi[11] = 0.044674560856694d
  164.   Ri[12] = 0.677872379632663d  &  Wi[12] = 0.047616658492490d
  165.   Ri[13] = 0.628867396776513d  &  Wi[13] = 0.050359035553854d
  166.   Ri[14] = 0.577224726083972d  &  Wi[14] = 0.052890189485193d
  167.   Ri[15] = 0.523160974722233d  &  Wi[15] = 0.055199503699984d
  168.   Ri[16] = 0.466902904750958d  &  Wi[16] = 0.057277292100403d
  169.   Ri[17] = 0.408686481990716d  &  Wi[17] = 0.059114839698395d
  170.   Ri[18] = 0.348755886292160d  &  Wi[18] = 0.060704439165893d
  171.   Ri[19] = 0.287362487355455d  &  Wi[19] = 0.062039423159892d
  172.   Ri[20] = 0.224763790394689d  &  Wi[20] = 0.063114192286254d
  173.   Ri[21] = 0.161222356068891d  &  Wi[21] = 0.063924238584648d
  174.   Ri[22] = 0.097004699209462d  &  Wi[22] = 0.064466164435950d
  175.   Ri[23] = 0.032380170962869d  &  Wi[23] = 0.064737696812683d
  176.   Ri[INDGEN(Pts/2) + (Pts/2)]  = -Ri[(Pts/2) - INDGEN(Pts/2) -1]
  177.   Wi[INDGEN(Pts/2) + (Pts/2)]  =  Wi[(Pts/2) - INDGEN(Pts/2) -1]
  178. endif else if Pts eq 96 then begin
  179.   Ri     = DBLARR(Pts)         &  Wi     = DBLARR(Pts)
  180.   Ri[0]  = 0.999689503883230d  &  Wi[0]  = 0.000796792065552d
  181.   Ri[1]  = 0.998364375863181d  &  Wi[1]  = 0.001853960788946d
  182.   Ri[2]  = 0.995981842987209d  &  Wi[2]  = 0.002910731817934d
  183.   Ri[3]  = 0.992543900323762d  &  Wi[3]  = 0.003964554338444d
  184.   Ri[4]  = 0.988054126329623d  &  Wi[4]  = 0.005014202742927d
  185.   Ri[5]  = 0.982517263563014d  &  Wi[5]  = 0.006058545504235d
  186.   Ri[6]  = 0.975939174585136d  &  Wi[6]  = 0.007096470791153d
  187.   Ri[7]  = 0.968326828463264d  &  Wi[7]  = 0.008126876925698d
  188.   Ri[8]  = 0.959688291448742d  &  Wi[8]  = 0.009148671230783d
  189.   Ri[9]  = 0.950032717784437d  &  Wi[9]  = 0.010160770535008d
  190.   Ri[10] = 0.939370339752755d  &  Wi[10] = 0.011162102099838d
  191.   Ri[11] = 0.927712456722308d  &  Wi[11] = 0.012151604671088d
  192.   Ri[12] = 0.915071423120898d  &  Wi[12] = 0.013128229566961d
  193.   Ri[13] = 0.901460635315852d  &  Wi[13] = 0.014090941772314d
  194.   Ri[14] = 0.886894517402420d  &  Wi[14] = 0.015038721026994d
  195.   Ri[15] = 0.871388505909296d  &  Wi[15] = 0.015970562902562d
  196.   Ri[16] = 0.854959033434601d  &  Wi[16] = 0.016885479864245d
  197.   Ri[17] = 0.837623511228187d  &  Wi[17] = 0.017782502316045d
  198.   Ri[18] = 0.819400310737931d  &  Wi[18] = 0.018660679627411d
  199.   Ri[19] = 0.800308744139140d  &  Wi[19] = 0.019519081140145d
  200.   Ri[20] = 0.780369043867433d  &  Wi[20] = 0.020356797154333d
  201.   Ri[21] = 0.759602341176647d  &  Wi[21] = 0.021172939892191d
  202.   Ri[22] = 0.738030643744400d  &  Wi[22] = 0.021966644438744d
  203.   Ri[23] = 0.715676812348967d  &  Wi[23] = 0.022737069658329d
  204.   Ri[24] = 0.692564536642171d  &  Wi[24] = 0.023483399085926d
  205.   Ri[25] = 0.668718310043916d  &  Wi[25] = 0.024204841792364d
  206.   Ri[26] = 0.644163403784967d  &  Wi[26] = 0.024900633222483d
  207.   Ri[27] = 0.618925840125468d  &  Wi[27] = 0.025570036005349d
  208.   Ri[28] = 0.593032364777572d  &  Wi[28] = 0.026212340735672d
  209.   Ri[29] = 0.566510418561397d  &  Wi[29] = 0.026826866725591d
  210.   Ri[30] = 0.539388108324357d  &  Wi[30] = 0.027412962726029d
  211.   Ri[31] = 0.511694177154667d  &  Wi[31] = 0.027970007616848d
  212.   Ri[32] = 0.483457973920596d  &  Wi[32] = 0.028497411065085d
  213.   Ri[33] = 0.454709422167743d  &  Wi[33] = 0.028994614150555d
  214.   Ri[34] = 0.425478988407300d  &  Wi[34] = 0.029461089958167d
  215.   Ri[35] = 0.395797649828908d  &  Wi[35] = 0.029896344136328d
  216.   Ri[36] = 0.365696861472313d  &  Wi[36] = 0.030299915420827d
  217.   Ri[37] = 0.335208522892625d  &  Wi[37] = 0.030671376123669d
  218.   Ri[38] = 0.304364944354496d  &  Wi[38] = 0.031010332586313d
  219.   Ri[39] = 0.273198812591049d  &  Wi[39] = 0.031316425596861d
  220.   Ri[40] = 0.241743156163840d  &  Wi[40] = 0.031589330770727d
  221.   Ri[41] = 0.210031310460567d  &  Wi[41] = 0.031828758894411d
  222.   Ri[42] = 0.178096882367618d  &  Wi[42] = 0.032034456231992d
  223.   Ri[43] = 0.145973714654896d  &  Wi[43] = 0.032206204794030d
  224.   Ri[44] = 0.113695850110665d  &  Wi[44] = 0.032343822568575d
  225.   Ri[45] = 0.081297495464425d  &  Wi[45] = 0.032447163714064d
  226.   Ri[46] = 0.048812985136049d  &  Wi[46] = 0.032516118713868d
  227.   Ri[47] = 0.016276744849602d  &  Wi[47] = 0.032550614492363d
  228.   Ri[INDGEN(Pts/2) + (Pts/2)]  = -Ri[(Pts/2) - INDGEN(Pts/2) -1]
  229.   Wi[INDGEN(Pts/2) + (Pts/2)]  =  Wi[(Pts/2) - INDGEN(Pts/2) -1]
  230. endif else MESSAGE, "Pts parameter must be 6, 10, 20, 48 or 96."
  231.  
  232.   TypeAB = SIZE(AB_Limits)
  233.  
  234.   ;If the DOUBLE keyword is not set then the internal precision and
  235.   ;result are identical to the type of input.
  236.   if N_ELEMENTS(Double) eq 0 then $
  237.     Double = (TypeAB[TypeAB[0]+1] eq 5) 
  238.  
  239.   if Double eq 0 then begin
  240.     Ri = FLOAT(Ri) & Wi = FLOAT(Wi)
  241.   endif
  242.  
  243.   nLimit = [AB_Limits[1] - AB_Limits[0], AB_Limits[1] + AB_Limits[0]] / 2.0
  244.   Aj = 0.0
  245.  
  246.   if KEYWORD_SET(Order) eq 0 then begin
  247.     ;Compute using a dy-dx order of integration.
  248.     for i = 0, Pts-1 do begin
  249.       X  = nLimit[0] * Ri[i] + nLimit[1]
  250.       fLimit = CALL_FUNCTION(PQ_Limits, X)
  251.       kF = (fLimit[1] - fLimit[0]) / 2.0
  252.       kS = (fLimit[1] + fLimit[0]) / 2.0
  253.       jX = TOTAL(Wi * CALL_FUNCTION(Fxy, X, kF*Ri+kS), Double = Double)
  254.       Aj = Aj + (Wi[i] * kF * jX)
  255.     endfor
  256.   endif else if Order eq 1 then begin 
  257.     ;Compute using a dx-dy order of integration.
  258.     for i = 0, Pts-1 do begin
  259.       Y  = nLimit[0] * Ri[i] + nLimit[1]
  260.       fLimit = CALL_FUNCTION(PQ_Limits, Y)
  261.       kF = (fLimit[1] - fLimit[0]) / 2.0
  262.       kS = (fLimit[1] + fLimit[0]) / 2.0
  263.       jY = TOTAL(Wi * CALL_FUNCTION(Fxy, kF*Ri+kS, Y), Double = Double)
  264.       Aj = Aj + (Wi[i] * kF * jY)
  265.     endfor
  266.   endif else MESSAGE, "Order keyword must be 0 or 1."
  267.  
  268.   if Double eq 0 then RETURN, FLOAT(Aj * nLimit[0]) else $
  269.     RETURN, (Aj * nLimit[0])
  270.  
  271. end
  272.